Thermalization in SU(3) gauge theory after a deconfining quench 
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We determine the time evolution of fluctuations of the Polyakov loop after a quench into the 
deconfined phase of SU(3) gauge theory from a simple classical relativistic Lagrangian. We com- 
pare the structure factors, which indicate spinodal decomposition followed by relaxation, to those 
obtained via Markov Chain Monte Carlo techniques in SU(3) lattice gauge theory. We find that 
the time when the structure factor peaks diverges like ~ in the long- wavelength limit. This 

is due to formation of competing Z(3) domains for configurations where the Polyakov loop exhibits 
non-perturbatively large variations in space, which delay thermalization of long wavelength modes. 
For realistic temperatures, and away from the extreme weak-coupling limit, we find that even modes 
with k on the order of T experience delayed thermalization. Relaxation times of very long wave- 
length modes are found to be on the order of the size of the system; thus, the dynamics of competing 
domains should accompany the hydrodynamic description of the deconfined vacuum. 



I. INTRODUCTION 

Relativistic Heavy-Ion Collision (RHIC) experiments 
carried out at Brookhaven National Laboratory (BNL) 
provide support for the existence of a quark-gluon plasma 
(QGP) phase of QCD In this phase color charges 
are liberated (deconfined), contrary to the low temper- 
ature phase that confines color charges inside colorless 
objects such as hadrons or glueballs. The existence of 
confined and deconfined phases has been demonstrated 
numerically in Lattice Gauge Theory (LOT) studies 0. 
Lattice QCD predicts a first-order phase transition for 
SU(3) pure gauge theory Q that turns into a crossover 
for full QCD [3| with physical quark masses. 

A collision of two heavy nuclei at high energy releases 
a large number of gluons from the wave functions of the 
colliding nuclei f^. Those gluons interact and eventu- 
ally form a thermalized QCD plasma with a temperature 
in excess of the critical temperature for deconfinement. 
Complete (and consistent) theoretical understanding of 
the thermalization process is presently lacking. Baier, 
Mueller, Schiff and Son developed the so-called "bottom- 
up" approach which is a framework for understand- 
ing the processes leading to thermalization and for cal- 
culating the thermalization time of the QCD medium 
as well as its initial temperature (see, however, the cri- 
tique in ref. Q)- The "bottom- up" approach is based 
on solving the Boltzmann equation for quasi-particles in 
a trivial vacuum and neglects the structure of the de- 
confined phase of the non-Abclian gauge theory arising 
from the Z(3) center symmetry discovered by t' Hooft 
and Polyakov [8]. Below, we shall show that in a model 
which allows for non-perturbatively large variations of 
Polyakov loops in space, domain walls form which sep- 



arate regions of different Z(3) orientation The com- 
petition between such domains affects the thermalization 
of long- wavelength modes of the Polyakov loop. We shall 
also provide some model estimates for the relevant wave 
lengths and time scales. 

We adopt a simplified picture where a relativistic 
heavy-ion collision is viewed as a quench that instantly 
heats the system to a temperature above the deconfining 
temperature. The response of Polyakov loop Structure 
Factors (SFs) to such a heating quench has been studied 
in SU(3) LGT by Markov Chain Monte Carlo (MCMC) 
simulations for Glauber (dissipative) dynamics [lO| . As 
unambiguous signals for the transition one finds a dy- 
namical growth of the SF, reaching maxima which scale 
approximately with the volume, a behavior often charac- 
terized as spinodal decomposition . 

Glauber (model A) dynamics 12] imitates the ther- 
mal fluctuation of Nature. It is therefore expected to 
describe the dissipative features of the transition from 
one equilibrium ensemble to another well. MCMC up- 
dating with a Metropolis or heatbath algorithm falls into 
the universality class of model A. Such SU(3) MCMC 
simulations converge to 3D equilibrium ensembles, which 
are the same as in Minkowski space, because the fourth 
extension of the Euclidean lattice serves only to define 
the temperature. The major drawback of Glauber dy- 
namics is not the 4D Euclidean formulation but that it 
is non-relativistic and, more importantly, that one does 
not know how to connect the MCMC updating step to a 
physical time scale. 

It was stressed in [13] that relaxation to the vacuum 
ensemble at high temperature becomes feasible only af- 
ter each SF has overcome its maximum value. For 
SU(3) gauge theory this relaxation time diverges with 
increasing system size due to competing order-order do- 
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mains with different Z(3) center group triahty, which 
are similar to order-order domains in a 3-state Potts 
model This divergence is well known in condensed 
matter physics 'l5'|. Hence, one should not a-priori ex- 
clude the possibility that, under heating, the long wave- 
length modes in the system do not equilibrate but instead 
get stuck in the neighborhood of the SF maxima. 

In this paper we explore heating quenches into the 
deconfined phase within Pisarski's effective model of 
Polyakov loops [l^ (see also ref. for similar effec- 
tive potential models), which we compare to those from 
Glauber dynamics. The solutions of the effective theory 
are obtained through Molecular Dynamics (MD) simula- 
tions of hyperbolic, relativistic Minkowski dynamics 
for which the time scale of the updating step is known 
(simulations based on a Langevin approach with statis- 
tical noise and friction have also been performed, see 
ref. [H). 

In Sec. |TT] we summarize the model and explain the 
simulations. Sec. lIIII presents our numerical results. Con- 
clusions and outlook follow in Sec. IIVI 

II. EFFECTIVE MODEL OF POLYAKOV LOOPS 

In Pisarski's model fl^ the deconfined phase of a 
pure gauge theory is described as a condensate of 
Polyakov loops. The Z(3)-symmetric effective potential 
for Polyakov loops t (complex for SU(3)) with cubic and 
quartic interactions takes the form 

m = - 1(^^ + {tf) + h ■ 

(1) 

The energy scale is set by T"*, and the mass coefficient 
b2 — b2{T) is temperature dependent, while 63 and 64 
are constants. These couplings can be chosen so that 
they reproduce lattice data for the SU(3) pressure and 
energy density above Tc- A reasonable fit follows with 
64 = 0.61 r*, 63 = 2.0/r, and b2{T) = ((1 - l.ll/x)(l + 
0.265/a;)2(l-|-0.300/j:)3-0.487)/r^ where x = T/T^ and 
r = 2.23. A plot is shown in fig. 1 of ref. [20| . 

To complete the effective theory in Minkowski space- 
time, a kinetic term has to be added: 

£ = t2 {Zt \dt£\^ + Zs \dd\^) - Vi£) . (2) 

Here, we assume a Lorentz-invariant form, Zt — Zs, and 
take the coefficient Zs from that for spatial variations 
of SU(3) Wilson hues [3, N^/g'^ with = 3. 

Thus the dynamics of Polyakov loops is in an intermedi- 
ate regime between very weak (Zs S> 1) and very strong 
(Zs <C 1) coupling. 

We employ a simulation procedure similar to the one of 
ref. [3 but focus on heating quenches of the system and 
not on its subsequent coohng [2l|; also, we consider here 
a static, non-expanding metric. Polyakov loop fields are 
defined on the sites of a spatial cubic lattice of size 
with periodic boundary conditions. They are initialized 
in the confined phase at time t = 0. Then the tem- 
perature entering the effective Lagrangian ^ is set to 
I 



a value Tf > Tc above the deconfinement transition Tc, 
where the Z(3) center symmetry is spontaneously broken. 
At temperatures Tf > Tc the effective potential takes the 
shape shown in Fig. [TJ In the plane below we show equal 
height contours. 




FIG. 1: Effective potential for the Polyakov loop £ from eq. ([T| 
at Tf /Tc = 2.0, shifted by a constant. Contours of equal 
heights are shown in the {£r,£i) plane. 

After the quench to Tf the system evolves in time ac- 
cording to the Euler-Lagrange equations derived from the 
Lagrangian ([2]) and rearranges itself to a new equilib- 
rium ensemble with non-zero Polyakov loop, which sig- 
nals symmetry breaking in the infinite volume limit. 

In the following we describe how the initial field con- 
figurations are constructed. 

A. Initialization of the field 

Introducing real and imaginary parts, we write the 
Polyakov loop as 

£ = i^ + ii^, (3) 

We split the fields into a long- and a short wavelength 
part, 

+ 6£r , £^^ei + 5£i . (4) 

The system is initialized in the unstable confined phase, 
thus 

ir{x, t = 0) = Z,(f , < = 0) = . (5) 

The initial fluctuations are assumed to be Gaussian, 

P[£r]-e^p^~^^ , P[4]~exp(^-^) . (6) 

Consequently, {5£r) = {S£^) = {5£.r5£,) = 0, {5£l) = 
and {5£l) = af hold. 
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B. Counterterms in the equations of motion 



In terms of ^r and ^i, the potential can be written as: 



(7) 



Using eq. ^ and averaging over the distribution of initial fluctuations (O up to order 5t'^ gives the one-loop correction 
to the effective potential: 



1 



■2/r/,2\ 



+ const 



(8) 



The additional second term, which arises from the classical initial fluctuations, needs to be subtracted in order to 
restore the original potential ([7]) for the long wavelength modes of the Polyakov loop. The equations of motion 
including these counterterms are given by 



1 r)V h 



(9) 



where we replaced i with £ because the difference in the equations of motion is of order 61^ and can be neglected. In 
the simulation we take a'^ — of = so that the counterterms simplify to 



1 av 



(10) 
(11) 



The numerical results shown below were obtained with 
a = 0.04; a = 0.08 gives similar results (our statistical 
errors are larger than the systematic errors due to varia- 
tion of the initial fluctuations within this range). We do 
not have a good quantitative estimate for the magnitude 
of fluctuations, but it appears reasonable to expect that 
over a length scale of order 1/Tc (see the next section) 
they should be small compared to unity. On the other 
hand, the value a = 0.04 is still sufficiently large to al- 
low for a relatively rapid onset of the domain formation 
process, as will be seen in section Hill 



C. Coarse graining 

A physical length scale is introduced through coarse 
graining of the initial field configuration. For example, a 
simple algorithm amounts to replacing the field £{x) at a 
given site by an average over a subvolume (box) of size 

£{x) ^ l'{x) = ^ E ^(^"0 ■ (12) 

^ff S'Gbox 



In physical units, the correlation length at should 
be on the order of 1/Tc (away from the extreme weak- 
coupling limit), and we set the lattice spacing a by 
aN^g = l/Tc. Dividing eqs. (jlOllip by T;^ shows that 
a different scale for the initial correlation length corre- 
sponds to rescaling Zg. 

The coarse-graining procedure does not affect the long 
wavelength part of i{x) but reduces the fluctuations, 

{Se'^) < . (13) 

Therefore, the counterterms in (jlOllip no longer match 
the state of the system after coarse-graining. One needs 
to restore the desired fluctuations by rescaling the initial 
flelds to 

6e"{x) = se'{x) , , (14) 

so that {51""^) = 0-2. These flelds are taken as initial 
configuration which is then propagated in time by solving 
the equations of motion. 



4 



D. Dynamics after a temperature quench 

To perform a quench the temperature is set to a vahie 
Tf in the deconfined phase. Then we use the leapfrog 
algorithm [13] to integrate the Euler-Lagrange equations 
(fTO|) . (fTTj) in time with the initial conditions described 
previously. At time the structure function is defined 
by the Fourier transformation of the Polyakov loop: 



F{k,t) 



~i k X 



£ix,t] 



(15) 



For a fixed value of k, F{k,t) is called SF. SFs are our 
primary observables. In what follows, we label SFs F„(t) 
similarly as in [l^: 



k = n2Tr/Ls , n = 1 
n = 2 
n — 3 



n 
n 
n 



(1,0,0), n' 
(1,1,0), 



1, (16) 

2, (17) 
(1,1,1), n^^3. (18) 




200 400 600 

t [monte carlo] 



FIG. 2: Structure factors for Glauber dynamics [lOj, t is 
Monte Carlo time in sweeps (4 X 64^ lattice, Tf /Tc = 1.57, 
lattice size Ls = 12.1 fm). Fn corresponds to the n*'^ lattice 
mode with = 2n^/Ls = V^x 102.1 MeV. 



where Ls = aNg denotes the size of the lattice in 
physical units. Note the relation = ^-K^fnjLs for 
n = 1, 2, 3. Measurements for n = \ include the permu- 
tations (0, 1, 0), (0, 0, 1) and for n = 2 the permutations 
(1,0,1), (0,1,1). 

III. NUMERICAL RESULTS 

We quench to several temperatures in the deconfining 
phase Ty/Tc = 1.50, 1.75, 2.00, 2.25, 2.50 on lattices with 
spatial extent = 40, 48, 64, 80, 96. Periodic bound- 
ary conditions are applied and averaging is done over 
ensembles of 200 replica. Our length scale set by coarse- 
graining is aNcg — 1/Tc = 0.736 fm, corresponding to the 
SU(3) phase transition temperature of Tc = 260 MeV [1]. 
Using A^cg = 4 for the correlation length the lattice spac- 
ing 

a = 0.184fm (19) 

follows. Physical volumes L\ — {aNg)^ in our simu- 
lations are (7.4)3, (8 g)3^ (11.8)^, (14.7)3 and (17.7)^ 
fm^. To reduce finite size effects we take lattices that 
accommodate at least 10 correlation lengths Ncg. When 
we study different physical volumes, Ncg has to be the 
same for all lattices and Ncg = 4 is a reasonable value. 
With, say, Ncg = 6 we would have to work on larger 
TVs = 60, • • • , 144 lattices. 

In Fig. [2] we present several SF modes from our Glauber 
dynamics study fld\ on a 4 x 64^ lattice (quench to 
(3 = 5.92, corresponding to Tf — 1.57 Tc, average over 
170 replica) for comparison with SF modes from the 
Minkowski dynamics on a 64^ lattice (quench to Tf = 
1.50 Tc, average over 200 rephca) presented in Fig. O 
The scale on the vertical axis of Fig. [2] differs from that 
of Fig. [3] because the former has been determined from 
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139(10) 
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TABLE I: -Fi,max for different volumes and temperatures 
(Minkowski dynamics). 

the bare Polyakov loop while the effective Lagrangian ^ 
deals with the renormalized loop. The renormalization 
constant for the Polyakov loop [23] amounts to a con- 
stant multiplying the SF. Since we are not interested in 
this renormalization here, this constant is of minor im- 
portance. 

Qualitatively, the SFs display the same behavior: an 
initial exponential growth is followed by equilibration af- 
ter the lowest SF mode reaches its maximum. As for 
Glauber dynamics, we interpret this as formation of com- 
peting order-order domains between regions of different 
Z(3) triality [9], whose equilibration takes a long time. 
For the Minkowskian dynamics the maxima of the lowest 
SF mode Fi are compiled in Table IH In the normaliza- 
tion of eq. (|15p they scale with volume in the same way 
as for Glauber dynamics, provided the volumes are large 
enough. Fits to the form -Fi^max = + ^iL^ are shown 
in Fig. m Interestingly, both Figs. [2] and [3] show that 
the above-mentioned effects due to non-perturbatively 
large variations of the Z(3) phase in space are visible 
even for modes with k ^ T/2. This may be related to 
Z(3) domain walls forming just after the quench being 
quite broad (disordered phase) and domains not much 
larger than 1/T. 

Since the kinetic term in the Lagrangian Q is as- 
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10 15 20 25 30 35 
t [fm/c] 



FIG. 3: Structure factors for Minkowski dynamics, t is real 
time (64^ lattice, Tf /Tc = 1.5, lattice size Ls = 11.8 fm). F„ 
corresponds to the n**^ lattice mode with |fc| = 2Tv^/n/Ls = 
^/nx 105.3 MeV. 




FIG. 4: Fits of -Fi.max for Minkowski dynamics. 



sumed to be Lorentz invariant, units for time are the 
same as for length (apart from the speed of hght fac- 
tor c). When integrating hyperbohc equations, one uses 
time steps smaller than the spatial lattice spacing a. We 
chose At /a — 0.01 and, therefore, in physical units 



At = 0.00184 fm/c 



(20) 



We ran trajectories from 15000 to 25000 time steps cor- 
responding to a range from 27.6 fm/c to 46 fm/c. 

For the case shown in Fig. [3] the structure factor for 
the first mode takes about 8 fm/c to reach its maximum, 
and another « 20 fm/c until that mode equilibrates. On 
the other hand, the second and third modes grow for a 
shorter period of time and subsequently equilibrate more 
rapidly. Note that there is an initial lag of w 2.5 fm/c, 
where growth of SFs is only visible on a logarithmic scale; 
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TABLE II: Times fmax [fm/c] where -Fi(t) peaks, for different 
volumes and temperatures (Minkowski dynamics). 




FIG. 5: Individual trajectories after a quench to Tf = 2.0 Tc 
on a 64^ lattice (Minkowski dynamics). Dashed lines show 
equipotential levels from Fig. [1] The circle about the origin 
indicates the average tmax time. 



the precise time for the onset of growth may be sensitive 
to the spectrum and magnitude of initial fluctuations. 

In Table ini we compile the times tmax needed to reach 
the maxima of Fi for different volumes and quench tem- 
peratures. Evidently, the times before the lowest modes 
of the system equilibrate can be quite large, in c = 1 units 
on the order of the size of the system. Note also that imax 
increases with Tf /Tc because the barriers between order- 
order domains grow higher and are more difficult to over- 
come. This is expected to change for full QCD with light 
quarks where the Z(3) symmetry is broken explicitly. 

Typical trajectories of the volume-averaged Polyakov 
loop are shown in Fig. together with contours of the 
potential from Fig. [TJ Some of the configurations in the 
ensemble relax directly towards one of the Z(3) minima, 
while others get stuck in-between for rather long times. 
The latter situation corresponds to the case when there 
are competing domains of approximately equal size: the 
relaxation of the long-wavelength modes is delayed when 
the Polyakov loop exhibits non-perturbatively large vari- 
ations in space. The circle of radius ctmax about the ori- 
gin indicates the average time imax- Note that for some 
trajectories the volume-averaged Polyakov loop (the total 
"magnetization" ) can at imax still be far from one of the 
minima of the potential. Individual configurations ex- 
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FIG. 6: Distribution of fmax after a quench to Tf — 2.0 Tc 
on a 64^ lattice for an ensemble of 1000 replica (Minkowski 
dynamics) . 
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FIG. 7: Fits of tmax for Minkowski dynamics. 



hibit rather large fluctuations about the mean, as shown 
in Fig. [H 

For Glauber dynamics it is known that for suffi- 
ciently large systems and corresponding wavelengths imax 
scales as 

imax = ao + ai . (21) 

For Minkowski dynamics fits to this form are satisfactory, 
too, as shown in Fig. [T] Thus, in the long wavelength 
limit the equilibration time for a mode with wave vector fc 
is proportional to The behavior ((2T|) could be used 

to scale ^max for F\ from table HI] into the corresponding 
^max for higher modes via the relation fc„ = L^. 
For Tf/Tc = 2.0, for example, we obtain ai = 0.044 fm"^ 
(and Oq = 2.6 fm, which corresponds to the "lag" men- 
tioned above). However, given the rather large statistical 
error bars shown in Fig. [7] we cannot exclude a linear de- 



pendence, 

^inax = ao + ^1 -^s J (22) 

with or without an initial "lag" ao, either. 

IV. CONCLUSIONS AND OUTLOOK 

The particle spectrum of a quantum field theory at fi- 
nite temperature builds upon the phase of its vacuum en- 
semble, and the structure of the vacuum is fundamental 
for the dynamics of relaxation processes to equilibrium. 
In SU(N) gauge theory, non-perturbatively large varia- 
tions of the Z(N) phase within domain walls arise during 
the conversion from a confined to a deconfined vacuum 
ensemble. We find that the Minkowskian dynamics of a 
simple model for SU(3) Polyakov loops (which incorpo- 
rates the Z(3) structure of the deconfined vacuum) re- 
produces the qualitative features of a previous study [loj 
of Glauber dynamics within LGT reasonably well. Heat- 
ing above Tc drives the SU(3) gauge theory system from 
the disordered into the ordered phase. The initial pro- 
cess of spinodal decomposition is signaled by an exponen- 
tial growth of the SFs. Ordering of the system proceeds 
through formation of domains of different Z(3) triality 
until one of them eventually occupies the whole system 
(see Fig. [SI . 

For realistic, sub-asymptotic temperatures as rele- 
vant for high-energy heavy-ion collisions we find that 
the rather slow dynamics of competing domains (of the 
Polyakov loop) delays thermalization of modes with k 
nearly up to T. In future, it would be interesting to 
quantify the contribution of these modes to the bulk vis- 
cosity. 

Our model for Minkowskian dynamics allows us to es- 
timate a physical time scale for the vacuum conversion 
process. The relaxation times found in this way for SU(3) 
pure gauge theory increase as ~ (for sufficiently 

low k) and are estimated to be on the order of the size 
Ls of the system for k = in/Ls and Lg ~ 10 fm. The 
hydrodynamic equations describing the evolution of long- 
wavelength perturbations in a deconfined medium should 
therefore be extended to account for the dynamics of 
competing domains of the Polyakov loop. 

The model for the dynamics of Polyakov loops could 
be improved in many ways. As an example, ([1]) deals 
only with the trace of the Polyakov loop in the funda- 
mental representation and neglects magnetic fields. Ef- 
fective Lagrangians which include all degrees of freedom 
of SU(N) loops as well as magnetic fields (in the static 
limit) have been proposed recently [23|. Also, the ki- 
netic term from ^ for non-equilibrium configurations, 
for which wc have assumed a Lorcntz-invariant form, is 
in fact not known. Further, one may worry about the 
importance of quantum corrections to the dynamics of 
domain walls. Finally, it would be interesting to include 
the effects due to dynamical quarks, which break the cen- 
ter symmetry explicitly. They reduce the SF maxima and 
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decrease the relaxation time, as has been demonstrated 
quahtatively for the 3D 3-state Potts model in an 
external magnetic field, while counter effects may come 
from the accompanying decrease of the transition tem- 
perature Tc and the transformation of the phase transi- 
tion into a crossover [ij. 

To improve the simulations, one could drop the as- 
sumption of instantaneous quenching and replace it by 
the thermalization time of the hard modes. The spinodal 
decomposition process starts already during the quench 
when the central rapidity volume heats up, and it occurs 
in an expanding medium where the longitudinal wave 
vectors experience a red-shift ^]. Moreover, boundary 
conditions in heavy-ion collisions are not periodic as used 
in this paper; rather, the surrounding vacuum is confined 
and therefore disordered. This increases the width and 
the effective temperature for the deconfinement transi- 
tion in a volume-dependent way. In a SU(3) equilibrium 



study the increase in temperature was found to be in 
the range from 20% for a volume of (5fm)'^ to 5% for a 
volume of (10 fm)^ 

Due to all of the above caveats, our numbers should 
be viewed only as rough first estimates. Nevertheless, 
it appears that away from the extreme weak-coupling 
limit [Zs S> 1) the dynamics of competing domains 
will influence thermalization of long wavelength modes, 
perhaps even for k not very far below T, and hence 
cannot be neglected. 
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